An end-to-end workflow for multiplexed image processing and analysis
Compiled on 19 January 2023
Abstract
Highly multiplexed imaging has enabled the simultaneous spatial profiling of dozens of biological molecules in tissues at single-cell resolution. Extracting biologically relevant information such as the spatial distribution of cell phenotypes from multiplexed tissue imaging data involves a number of computational tasks, including image segmentation, feature extraction, and spatially-resolved single-cell analysis. Here, we present an end-to-end workflow for multiplexed tissue image processing and analysis, integrating a number of previously developed computational tools to enable these tasks in a user-friendly and customizable fashion. For data quality assessment, we highlight the utility of napari-imc for interactively inspecting raw imaging data and the cytomapper R/Bioconductor package for image visualization in R. Raw data preprocessing, image segmentation and feature extraction are performed using the steinbock toolkit. We showcase two alternative approaches for segmenting cells based on supervised pixel classification and pre-trained deep learning models. The extracted single-cell data is then read, processed and analyzed in R. The protocol describes the use of common data containers, facilitating the application of a number of R/Bioconductor packages for dimensionality reduction, single-cell visualization and phenotyping. We provide instructions to perform spatially-resolved single-cell analysis including community analysis, cellular neighborhood detection and cell-cell interaction testing using the imcRtools R/Bioconductor package. Overall, we provide the protocol for researchers with little bioinformatics training, and data analysis can be completed within 5-6 hours, depending on the segmentation approach. An extended version of the workflow can be accessed at https://bodenmillergroup.github.io/IMCDataAnalysis/.
1 Materials
1.1 Data
The example data needed to run the protocol can be downloaded as follows:
options(timeout = 10000)
dir.create("data/steinbock/raw", recursive = TRUE)## Warning in dir.create("data/steinbock/raw", recursive = TRUE):
## 'data/steinbock/raw' already existsdownload.file("https://zenodo.org/record/7412972/files/panel.csv",
"data/steinbock/panel.csv")
download.file("https://zenodo.org/record/5949116/files/Patient1.zip",
"data/steinbock/raw/Patient1.zip")
download.file("https://zenodo.org/record/5949116/files/Patient2.zip",
"data/steinbock/raw/Patient2.zip")
download.file("https://zenodo.org/record/5949116/files/Patient3.zip",
"data/steinbock/raw/Patient3.zip")
download.file("https://zenodo.org/record/5949116/files/Patient4.zip",
"data/steinbock/raw/Patient4.zip")
download.file("https://zenodo.org/record/5949116/files/compensation.zip",
"data/compensation.zip")
unzip("data/compensation.zip", exdir="data", overwrite=TRUE)
unlink("data/compensation.zip")
download.file("https://zenodo.org/record/5949116/files/sample_metadata.xlsx",
destfile = "data/sample_metadata.xlsx")
download.file("https://zenodo.org/record/7432486/files/gated_cells.zip",
"data/gated_cells.zip")
unzip("data/gated_cells.zip", exdir="data", overwrite=TRUE)
unlink("data/gated_cells.zip")2 Equipment setup
2.1 Installation instructions
To install all needed R packages, run the following code:
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github(c("BodenmillerGroup/imcRtools",
"BodenmillerGroup/cytomapper"))
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("pheatmap", "viridis",
"tiff", "distill", "openxlsx", "ggrepel", "patchwork",
"mclust", "RColorBrewer", "uwot", "Rtsne", "caret",
"randomForest", "ggridges", "gridGraphics", "scales",
"CATALYST", "scuttle", "scater", "dittoSeq",
"tidyverse", "batchelor", "bluster","scran"))3 Procedure
3.1 Setting the steinbock alias
TODO
alias steinbock="docker run -v /path/to/IMCDataAnalysis/data/steinbock:/data -u $(id -u):$(id -g) ghcr.io/bodenmillergroup/steinbock:0.15.0"3.2 Preprocessing
In this protocol, data pre-processing refers to the extraction of multi-channel images from raw imaging data, and to preparing them for downstream processing. The required steps are dependent on the imaging technology; here, we showcase the pre-processing of raw IMC data which includes a hot pixel filtering step.
steinbock preprocess imc images --hpf 50## bash: no job control in this shell
## 2023-01-19 11:04:53,546 INFO steinbock - img/Patient4_005.tiff
## 2023-01-19 11:04:54,406 INFO steinbock - img/Patient4_006.tiff
## 2023-01-19 11:04:55,216 INFO steinbock - img/Patient4_007.tiff
## 2023-01-19 11:04:55,864 INFO steinbock - img/Patient4_008.tiff
## 2023-01-19 11:05:00,603 INFO steinbock - img/Patient3_001.tiff
## 2023-01-19 11:05:01,251 INFO steinbock - img/Patient3_002.tiff
## 2023-01-19 11:05:01,966 INFO steinbock - img/Patient3_003.tiff
## 2023-01-19 11:05:07,272 INFO steinbock - img/Patient2_001.tiff
## 2023-01-19 11:05:08,039 INFO steinbock - img/Patient2_002.tiff
## 2023-01-19 11:05:08,827 INFO steinbock - img/Patient2_003.tiff
## 2023-01-19 11:05:09,631 INFO steinbock - img/Patient2_004.tiff
## 2023-01-19 11:05:13,329 INFO steinbock - img/Patient1_001.tiff
## 2023-01-19 11:05:14,069 INFO steinbock - img/Patient1_002.tiff
## 2023-01-19 11:05:14,882 INFO steinbock - img/Patient1_003.tiff
## 2023-01-19 11:05:14,955 INFO steinbock - images.csvThe step took 0.57 minutes.
3.3 Image segmentation
Perform automatic deep learning-enabled single-cell segmentation using the pre-trained Mesmer neural network implemented in DeepCell. In the following command, channels will be min-max-normalized and mean-aggregated according to the deepcell column in the panel file.
steinbock segment deepcell --minmax## bash: no job control in this shell
## /opt/venv/lib/python3.8/site-packages/deepcell_toolbox/deep_watershed.py:179: FutureWarning: `selem` is a deprecated argument name for `h_maxima`. It will be removed in version 1.0. Please use `footprint` instead.
## markers = h_maxima(image=maxima,
## 2023-01-19 11:05:46,458 INFO steinbock - masks/Patient1_001.tiff
## 2023-01-19 11:05:59,741 INFO steinbock - masks/Patient1_002.tiff
## 2023-01-19 11:06:12,880 INFO steinbock - masks/Patient1_003.tiff
## 2023-01-19 11:06:25,566 INFO steinbock - masks/Patient2_001.tiff
## 2023-01-19 11:06:37,911 INFO steinbock - masks/Patient2_002.tiff
## 2023-01-19 11:06:51,265 INFO steinbock - masks/Patient2_003.tiff
## 2023-01-19 11:07:03,837 INFO steinbock - masks/Patient2_004.tiff
## 2023-01-19 11:07:17,476 INFO steinbock - masks/Patient3_001.tiff
## 2023-01-19 11:07:29,921 INFO steinbock - masks/Patient3_002.tiff
## 2023-01-19 11:07:43,243 INFO steinbock - masks/Patient3_003.tiff
## 2023-01-19 11:07:56,739 INFO steinbock - masks/Patient4_005.tiff
## 2023-01-19 11:08:10,194 INFO steinbock - masks/Patient4_006.tiff
## 2023-01-19 11:08:21,969 INFO steinbock - masks/Patient4_007.tiff
## 2023-01-19 11:08:33,911 INFO steinbock - masks/Patient4_008.tiffThe step took 3.34 minutes.
3.4 Single-cell data extraction
For each image, extract the mean pixel intensity per cell and marker. The resulting cell-level intensity values are stored as separate CSV files (one file per image):
steinbock measure intensities## bash: no job control in this shell
## 2023-01-19 11:08:38,573 INFO steinbock - intensities/Patient1_001.csv
## 2023-01-19 11:08:39,535 INFO steinbock - intensities/Patient1_002.csv
## 2023-01-19 11:08:40,633 INFO steinbock - intensities/Patient1_003.csv
## 2023-01-19 11:08:41,717 INFO steinbock - intensities/Patient2_001.csv
## 2023-01-19 11:08:42,675 INFO steinbock - intensities/Patient2_002.csv
## 2023-01-19 11:08:43,554 INFO steinbock - intensities/Patient2_003.csv
## 2023-01-19 11:08:44,705 INFO steinbock - intensities/Patient2_004.csv
## 2023-01-19 11:08:45,819 INFO steinbock - intensities/Patient3_001.csv
## 2023-01-19 11:08:47,775 INFO steinbock - intensities/Patient3_002.csv
## 2023-01-19 11:08:48,819 INFO steinbock - intensities/Patient3_003.csv
## 2023-01-19 11:08:49,799 INFO steinbock - intensities/Patient4_005.csv
## 2023-01-19 11:08:50,795 INFO steinbock - intensities/Patient4_006.csv
## 2023-01-19 11:08:51,670 INFO steinbock - intensities/Patient4_007.csv
## 2023-01-19 11:08:53,039 INFO steinbock - intensities/Patient4_008.csvThe step took 0.29 minutes.
For each image, extract standard morphological features (e.g., area, eccentricity) per cell. The resulting cell-level features are stored as separate CSV files (one file per image):
steinbock measure regionprops## bash: no job control in this shell
## 2023-01-19 11:08:56,102 INFO steinbock - regionprops/Patient1_001.csv
## 2023-01-19 11:08:57,284 INFO steinbock - regionprops/Patient1_002.csv
## 2023-01-19 11:08:58,549 INFO steinbock - regionprops/Patient1_003.csv
## 2023-01-19 11:08:59,672 INFO steinbock - regionprops/Patient2_001.csv
## 2023-01-19 11:09:00,847 INFO steinbock - regionprops/Patient2_002.csv
## 2023-01-19 11:09:01,812 INFO steinbock - regionprops/Patient2_003.csv
## 2023-01-19 11:09:03,127 INFO steinbock - regionprops/Patient2_004.csv
## 2023-01-19 11:09:04,532 INFO steinbock - regionprops/Patient3_001.csv
## 2023-01-19 11:09:05,714 INFO steinbock - regionprops/Patient3_002.csv
## 2023-01-19 11:09:06,975 INFO steinbock - regionprops/Patient3_003.csv
## 2023-01-19 11:09:08,022 INFO steinbock - regionprops/Patient4_005.csv
## 2023-01-19 11:09:09,487 INFO steinbock - regionprops/Patient4_006.csv
## 2023-01-19 11:09:10,576 INFO steinbock - regionprops/Patient4_007.csv
## 2023-01-19 11:09:11,637 INFO steinbock - regionprops/Patient4_008.csvThe step took 0.31 minutes.
In each image, detect cells in close spatial proximity. The resulting spatial cell graphs are stored as separate directed edge lists in CSV format (one file per image):
steinbock measure neighbors --type expansion --dmax 4## bash: no job control in this shell
## 2023-01-19 11:09:16,799 INFO steinbock - neighbors/Patient1_001.csv
## 2023-01-19 11:09:19,852 INFO steinbock - neighbors/Patient1_002.csv
## 2023-01-19 11:09:23,404 INFO steinbock - neighbors/Patient1_003.csv
## 2023-01-19 11:09:26,336 INFO steinbock - neighbors/Patient2_001.csv
## 2023-01-19 11:09:29,224 INFO steinbock - neighbors/Patient2_002.csv
## 2023-01-19 11:09:31,628 INFO steinbock - neighbors/Patient2_003.csv
## 2023-01-19 11:09:35,282 INFO steinbock - neighbors/Patient2_004.csv
## 2023-01-19 11:09:38,903 INFO steinbock - neighbors/Patient3_001.csv
## 2023-01-19 11:09:42,087 INFO steinbock - neighbors/Patient3_002.csv
## 2023-01-19 11:09:45,597 INFO steinbock - neighbors/Patient3_003.csv
## 2023-01-19 11:09:48,152 INFO steinbock - neighbors/Patient4_005.csv
## 2023-01-19 11:09:52,139 INFO steinbock - neighbors/Patient4_006.csv
## 2023-01-19 11:09:55,011 INFO steinbock - neighbors/Patient4_007.csv
## 2023-01-19 11:09:57,745 INFO steinbock - neighbors/Patient4_008.csvThe step took 0.77 minutes.
3.5 Reading in data
Read in the spatially-resolved single-cell data into R using the imcRtools package. For the rest of the protocol we will continue with the steinbock generated data.
library(imcRtools)
spe <- read_steinbock("data/steinbock/")The step took 0.57 minutes.
After reading in the single-cell data, the SpatialExperiment object needs to be further processed.
First, the column names are set based on the image name and the cell identifier.
The patient identifier and the region of interest (ROI) identifier are saved in the object as well as the cancer type, which can be read in from the provided data/sample_metadata.xlsx file.
For easy access later on, the channels containing biological variation are selected.
Finally, the mean pixel intensities per channel and cell are arsinh-transformed using a cofactor of 1.
library(openxlsx)
library(tidyverse)
colnames(spe) <- paste0(spe$sample_id, "_", spe$ObjectNumber)
meta <- read.xlsx("data/sample_metadata.xlsx")
spe$patient_id <- as.vector(str_extract_all(spe$sample_id, "Patient[1-4]",
simplify = TRUE))
spe$ROI <- as.vector(str_extract_all(spe$sample_id, "00[1-8]",
simplify = TRUE))
spe$indication <- meta$Indication[match(spe$patient_id, meta$Sample.ID)]
rowData(spe)$use_channel <- !grepl("DNA|Histone", rownames(spe))
assay(spe, "exprs") <- asinh(counts(spe)/1)The step took 0.04 minutes.
Read in multi-channel images as a CytoImageList container using the cytomapper package.
library(cytomapper)
images <- loadImages("data/steinbock/img/")
channelNames(images) <- rownames(spe)The step took 0.26 minutes.
Read in segmentation masks as a CytoImageList container.
masks <- loadImages("data/steinbock/masks/", as.is = TRUE)## All files in the provided location will be read in.The step took 0 minutes.
For downstream visualization and analysis tasks, additional metadata needs to be added to the CytoImageList objects storing the multi-channel images and segmentation masks. Here, individual images, segmentation masks and entries in the SpatialExperiment object are matched via the sample_id entry.
patient_id <- str_extract_all(names(images), "Patient[1-4]", simplify = TRUE)
indication <- meta$Indication[match(patient_id, meta$Sample.ID)]
mcols(images) <- mcols(masks) <- DataFrame(sample_id = names(images),
patient_id = patient_id,
indication = indication)The step took 0 minutes.
3.6 Spillover correction
Low signal spillover between neighbouring channels occurs when using technologies that rely on mass cytometry. Spillover is defined as a small proportion of the signal of a neighbouring channel which can be detected in the primary channel. As spillover is linear to the signal of the neighbouring channel it can be correct for using a previously described compensation approach. This phenomenon is IMC specific and the steps of the following section can be skipped when working with data generated by other multiplexed imaging technologies.
Read in data from the spillover slide for channel-to-channel spillover correction. The experimental procedure to create and acquire a spillover slide can be seen in Supplementary Note 2. As recommended by the CATALYST R/Bioconductor package, the pixel intensities are arsinh-transformed using a cofactor of 5.
sce <- readSCEfromTXT("data/compensation/")## Spotted channels: Y89, In113, In115, Pr141, Nd142, Nd143, Nd144, Nd145, Nd146, Sm147, Nd148, Sm149, Nd150, Eu151, Sm152, Eu153, Sm154, Gd155, Gd156, Gd158, Tb159, Gd160, Dy161, Dy162, Dy163, Dy164, Ho165, Er166, Er167, Er168, Tm169, Er170, Yb171, Yb172, Yb173, Yb174, Lu175, Yb176
## Acquired channels: Ar80, Y89, In113, In115, Xe131, Xe134, Ba136, La138, Pr141, Nd142, Nd143, Nd144, Nd145, Nd146, Sm147, Nd148, Sm149, Nd150, Eu151, Sm152, Eu153, Sm154, Gd155, Gd156, Gd158, Tb159, Gd160, Dy161, Dy162, Dy163, Dy164, Ho165, Er166, Er167, Er168, Tm169, Er170, Yb171, Yb172, Yb173, Yb174, Lu175, Yb176, Ir191, Ir193, Pt196, Pb206
## Channels spotted but not acquired:
## Channels acquired but not spotted: Ar80, Xe131, Xe134, Ba136, La138, Ir191, Ir193, Pt196, Pb206assay(sce, "exprs") <- asinh(counts(sce)/5)The step took 0.11 minutes.
Perform quality assessment of the spillover data by visualizing the median pixel intensity per channel and spotted metal.
plotSpotHeatmap(sce)
(optional) Perform pixel binning to increase median pixel intensity. This is only needed if pixel intensities are too low (median below ~200 counts).
sce2 <- binAcrossPixels(sce, bin_size = 10)The step took 0.2 minutes.
Filter incorrectly assigned pixels. The following step uses functions provided by the CATALYST package to “de-barcode” the pixels. Based on the intensity distribution of all channels, pixels are assigned to their corresponding barcode; here, this is the already known metal spot. This procedure identifies pixels that cannot be robustly assigned to the spotted metal. Pixels of such kind can be regarded as “noisy”, “background”, or “artifacts” that should be removed prior to spillover estimation. The spotted channels (bc_key) need to be specified. The general workflow for pixel de-barcoding is as follows:
- assign a preliminary metal mass to each pixel
- for each pixel, estimate a cutoff parameter for the distance between positive and negative pixel sets
- apply the estimated cutoffs to identify truly positive pixels
In cases where incorrect assignments occurred or where few pixels were measured for some spots, the imcRtools package exports a helper function to exclude pixels.
library(CATALYST)
bc_key <- as.numeric(unique(sce$sample_mass))
bc_key <- bc_key[order(bc_key)]
sce <- assignPrelim(sce, bc_key = bc_key)## Debarcoding data...## o ordering## o classifying events## Normalizing...## Computing deltas...sce <- estCutoffs(sce)
sce <- applyCutoffs(sce)
sce <- filterPixels(sce, minevents = 40, correct_pixels = TRUE)The step took 0.14 minutes.
Compute and store the spillover matrix using the CATALYST package.
sce <- computeSpillmat(sce)
sm <- metadata(sce)$spillover_matrixThe step took 0.03 minutes.
Perform single-cell data compensation using the CATALYST package. The compCytof function corrects channel-to-channel spillover directly on the single-cell intensities using the previously estimated spillover matrix. The isotope_list variable needs to be extended by isotopes that are not contained in this list provided by the CATALYST package. Visualization of marker intensities of neighboring channels (e.g., Yb173 and Yb174) before and after correction can be used to assess the spillover correction efficacy.
library(dittoSeq)
library(patchwork)
rowData(spe)$channel_name <- paste0(rowData(spe)$channel, "Di")
isotope_list <- CATALYST::isotope_list
isotope_list$Ar <- 80
spe <- compCytof(spe, sm,
transform = TRUE, cofactor = 1,
isotope_list = isotope_list,
overwrite = FALSE)
before <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
assay.x = "exprs", assay.y = "exprs") +
ggtitle("Before compensation")
after <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
assay.x = "compexprs", assay.y = "compexprs") +
ggtitle("After compensation")
before + after
assay(spe, "counts") <- assay(spe, "compcounts")
assay(spe, "exprs") <- assay(spe, "compexprs")
assay(spe, "compcounts") <- assay(spe, "compexprs") <- NULLThe step took 0.2 minutes.
Perform channel-to-channel spillover correction on multi-channel images. To this end, the previously computed spillover matrix needs to be adjusted to only retain channels that are stored in the multi-channel images. By visualizing neighboring channels, spillover correction efficacy can be assessed.
channelNames(images) <- rowData(spe)$channel_name
adapted_sm <- adaptSpillmat(sm, paste0(rowData(spe)$channel, "Di"),
isotope_list = isotope_list)## Compensation is likely to be inaccurate.
## Spill values for the following interactions
## have not been estimated:## Ir191Di -> Ir193Di## Ir193Di -> Ir191Diimages_comp <- compImage(images, adapted_sm)
plotPixels(images[5], colour_by = "Yb173Di",
image_title = list(text = "Yb173 (Ecad) - before",
position = "topleft"),
legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
plotPixels(images[5], colour_by = "Yb174Di",
image_title = list(text = "Yb174 (CD303) - before",
position = "topleft"),
legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
plotPixels(images_comp[5], colour_by = "Yb173Di",
image_title = list(text = "Yb173 (Ecad) - after",
position = "topleft"),
legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
plotPixels(images_comp[5], colour_by = "Yb174Di",
image_title = list(text = "Yb174 (CD303) - after",
position = "topleft"),
legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
channelNames(images_comp) <- rownames(spe)The step took 9.59 minutes.
3.7 Quality control
Outline cells on composite images for visual assessment of segmentation quality. For visualization purposes, we subset 3 images and outline all cells on composite images after channel normalization.
set.seed(20220118)
img_ids <- sample(seq_len(length(images_comp)), 3)
cur_images <- images_comp[img_ids]
cur_images <- normalize(cur_images, separateImages = TRUE)
cur_images <- normalize(cur_images, inputRange = c(0, 0.2))
plotPixels(cur_images,
mask = masks[img_ids],
img_id = "sample_id",
missing_colour = "white",
colour_by = c("CD163", "CD20", "CD3", "Ecad", "DNA1"),
colour = list(CD163 = c("black", "yellow"),
CD20 = c("black", "red"),
CD3 = c("black", "green"),
Ecad = c("black", "cyan"),
DNA1 = c("black", "blue")),
image_title = NULL,
legend = list(colour_by.title.cex = 0.9,
colour_by.labels.cex = 0.9))
Visualize the distribution of the cell area and filter out small cells.
colData(spe) %>%
as.data.frame() %>%
group_by(sample_id) %>%
ggplot() +
geom_boxplot(aes(sample_id, area)) +
theme_minimal(base_size = 15) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 15)) +
ylab("Cell area") + xlab("")
spe <- spe[,spe$area >= 5]Visualize the cell density per image.
colData(spe) %>%
as.data.frame() %>%
group_by(sample_id) %>%
summarize(cell_area = sum(area),
no_pixels = mean(width_px) * mean(height_px)) %>%
mutate(covered_area = cell_area / no_pixels) %>%
ggplot() +
geom_point(aes(sample_id, covered_area)) +
theme_minimal(base_size = 15) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 15)) +
ylim(c(0, 1)) +
ylab("% covered area") + xlab("")
Visualize staining differences between samples for selected markers. Together with the visualization of cells in low dimensions, this figure indicates sample-to-sample differences in marker expression.
multi_dittoPlot(spe, vars = c("HLADR", "CD3", "Ecad", "PDGFRb"),
group.by = "patient_id", plots = c("ridgeplot"),
assay = "exprs")## Picking joint bandwidth of 0.162## Picking joint bandwidth of 0.0982## Picking joint bandwidth of 0.151## Picking joint bandwidth of 0.089
Visualize low-dimensional embeddings of single cells. Here, we use the scater package to compute a Uniform Manifold Approximation and Projection (UMAP) embedding and visualize cells in low-dimensional space.
library(scater)## Loading required package: scuttleset.seed(220225)
spe <- runUMAP(spe, subset_row = rowData(spe)$use_channel,
exprs_values = "exprs")
dittoDimPlot(spe, var = "patient_id",
reduction.use = "UMAP", size = 0.2) +
ggtitle("Patient ID on UMAP")
The step took 0.58 minutes.
Perform batch correction to remove sample-to-sample differences. Here, we use the fastMNN method of the batchelor package.
library(batchelor)
set.seed(220228)
out <- fastMNN(spe, batch = spe$patient_id,
auto.merge = TRUE,
subset.row = rowData(spe)$use_channel,
assay.type = "exprs")
reducedDim(spe, "fastMNN") <- reducedDim(out, "corrected")
spe <- runUMAP(spe, dimred= "fastMNN", name = "UMAP_mnnCorrected")
dittoDimPlot(spe, var = "patient_id",
reduction.use = "UMAP_mnnCorrected", size = 0.2) +
ggtitle("Patient ID on UMAP after correction")
The step took 2.1 minutes.
3.8 Cell phenotyping
Define cellular phenotypes. For this, single cells can be clustered (A) or cells can be labelled via classification (B).
3.8.1 (A) Cell phenotyping via clustering.
Graph-based clustering is performed using functions from the bluster and scran R/Bioconductor packages. Alternatively, other approaches such as phenograph or FlowSOM can be used to cluster single cells.
Estimate optimal clustering parameters for graph-based clustering on the integrated cells after batch correction. We perform a sweep across possible combinations of clustering parameters, including the number of nearest neighbors to consider (k) and the edge weighting method (type). We keep the community detection algorithm (cluster.fun) fixed, as the Louvain method is one of the most commonly used algorithms for graph-based clustering. To assess cluster stability, we compute the mean silhouette width across all cells and select the cluster parameter combination with highest mean silhouette width.
library(bluster)
library(BiocParallel)
mat <- reducedDim(spe, "fastMNN")
set.seed(220825)
combinations <- clusterSweep(mat, BLUSPARAM=SNNGraphParam(),
k=c(10L, 20L),
type = c("rank", "jaccard"),
cluster.fun = "louvain")
sil <- vapply(as.list(combinations$clusters),
function(x) mean(approxSilhouette(mat, x)$width), 0)
ggplot(data.frame(method = names(sil),
sil = sil)) +
geom_point(aes(method, sil), size = 3) +
theme_classic(base_size = 15) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
xlab("Cluster parameter combination") +
ylab("Average silhouette width")
The step took 2.99 minutes.
Based on the selected parameters, cells are clustered using a graph-based algorithm. As observed above, the parameter setting for k=20 and type=”rank” should result in well-separated clusters. The cluster identifiers are then saved in the SpatialExperiment object.
library(scran)
clusters <- clusterCells(spe[rowData(spe)$use_channel,],
use.dimred = "fastMNN",
BLUSPARAM = SNNGraphParam(k=20,
cluster.fun = "louvain",
type = "rank"))
spe$nn_clusters <- clustersThe step took 1.32 minutes.
3.8.2 (B) Classification-based cell phenotyping
Gate individual cell types based on their marker expression. For this, the cytomapper package provides the cytomapperShiny function. Per image, cells are gated based on their marker expression in a hierarchical fashion to define expected cell types. The gated cells are then visualized as outlines on pseudo-coloured composite images. Once the correct cells are labeled, they can be downloaded as a SpatialExperiment object storing only the selected cells. During download, the cell label can be specified, which is stored in the cytomapper_CellLabel entry of the colData slot for later use in training a classifier.
if (interactive()){
cytomapperShiny(object = spe, mask = masks, image = images_comp,
cell_id = "ObjectNumber", img_id = "sample_id")
}Read in files containing the gated cells and concatenate them into a single SpatialExperiment object.
library(SingleCellExperiment)
label_files <- list.files("data/gated_cells",
full.names = TRUE, pattern = ".rds$")
spes <- lapply(label_files, readRDS)
concat_spe <- do.call("cbind", spes)Remove doublets and reassign tumor cells. As a result of the multi-step labeling approach, some cells may have been labeled multiple times. In cases where one cell was labeled as both tumor and immune cell, we keep the immune cell label, as these cells are most likely immune cells residing within the tumor. All other cells that were labeled multiple times are removed. Finally, the labels are stored in the main SpatialExperiment object.
cur_tab <- unclass(table(colnames(concat_spe),
concat_spe$cytomapper_CellLabel))
cur_labels <- rep("doublets", nrow(cur_tab))
names(cur_labels) <- rownames(cur_tab)
# Single assignments
single_index <- rowSums(cur_tab) == 1
cur_labels[single_index] <- colnames(cur_tab)[
apply(cur_tab[single_index,], 1, which.max)]
# Double assignment within the tumor
double_index <- rowSums(cur_tab) == 2 & cur_tab[,"Tumor"] == 1
no_tumor <- cur_tab[,colnames(cur_tab) != "Tumor"]
cur_labels[double_index] <- colnames(no_tumor)[
apply(no_tumor[double_index,], 1, which.max)]
# Remove doublets
cur_labels <- cur_labels[cur_labels != "doublets"]
# Transfer labels to SPE object
spe_labels <- rep("unlabeled", ncol(spe))
names(spe_labels) <- colnames(spe)
spe_labels[names(cur_labels)] <- cur_labels
spe$cell_labels <- spe_labelsTrain a random forest classifier for cell type classification of unlabelled cells. The cells are first split into labelled and unlabelled cells. We then perform a 75/25 split of the labelled cells to select training and testing datasets, respectively. Based on the training dataset, we perform a 5-fold cross validation to tune the random forest model parameter.
library(caret)
# Split between labeled and unlabeled cells
lab_spe <- spe[,spe$cell_labels != "unlabeled"]
unlab_spe <- spe[,spe$cell_labels == "unlabeled"]
# Randomly split into train and test data
set.seed(220925)
trainIndex <- createDataPartition(factor(lab_spe$cell_labels), p = 0.75)
train_spe <- lab_spe[,trainIndex$Resample1]
test_spe <- lab_spe[,-trainIndex$Resample1]
fitControl <- trainControl(method = "cv",
number = 5)
cur_mat <- t(assay(train_spe, "exprs")[rowData(train_spe)$use_channel,])
rffit <- train(x = cur_mat,
y = factor(train_spe$cell_labels),
method = "rf", ntree = 1000,
tuneLength = 5,
trControl = fitControl)The step took 10.42 minutes.
Assess the classifier performance by computing the confusion matrix of the test dataset. The confusionMatrix function compares the predicted cell labels to the ground truth cell labels and computes a number of performance metrics. A high sensitivity and a high specificity for each cell type label is to be desired.
cur_mat <- t(assay(test_spe, "exprs")[rowData(test_spe)$use_channel,])
cur_pred <- predict(rffit,
newdata = cur_mat)
confusionMatrix(data = cur_pred,
reference = factor(test_spe$cell_labels),
mode = "everything")## Confusion Matrix and Statistics
##
## Reference
## Prediction Bcell BnTcell CD4 CD8 Myeloid Neutrophil Plasma_cell Stroma
## Bcell 188 0 0 0 0 0 7 0
## BnTcell 1 423 0 0 0 0 0 0
## CD4 0 0 164 0 0 0 7 0
## CD8 0 0 0 200 0 0 7 0
## Myeloid 0 0 4 2 437 0 0 0
## Neutrophil 0 0 0 0 0 32 0 0
## Plasma_cell 1 0 2 2 0 0 155 0
## Stroma 0 0 0 0 0 0 0 110
## Treg 0 0 0 0 0 0 2 0
## Tumor 5 2 2 1 0 1 1 0
## Reference
## Prediction Treg Tumor
## Bcell 0 1
## BnTcell 0 0
## CD4 1 2
## CD8 0 6
## Myeloid 0 0
## Neutrophil 0 2
## Plasma_cell 1 1
## Stroma 0 0
## Treg 87 2
## Tumor 1 1485
##
## Overall Statistics
##
## Accuracy : 0.9809
## 95% CI : (0.9756, 0.9852)
## No Information Rate : 0.4481
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9745
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Bcell Class: BnTcell Class: CD4 Class: CD8
## Sensitivity 0.96410 0.9953 0.95349 0.97561
## Specificity 0.99746 0.9997 0.99685 0.99586
## Pos Pred Value 0.95918 0.9976 0.94253 0.93897
## Neg Pred Value 0.99778 0.9993 0.99748 0.99840
## Precision 0.95918 0.9976 0.94253 0.93897
## Recall 0.96410 0.9953 0.95349 0.97561
## F1 0.96164 0.9965 0.94798 0.95694
## Prevalence 0.05830 0.1271 0.05142 0.06129
## Detection Rate 0.05620 0.1265 0.04903 0.05979
## Detection Prevalence 0.05859 0.1268 0.05202 0.06368
## Balanced Accuracy 0.98078 0.9975 0.97517 0.98573
## Class: Myeloid Class: Neutrophil Class: Plasma_cell
## Sensitivity 1.0000 0.969697 0.86592
## Specificity 0.9979 0.999396 0.99779
## Pos Pred Value 0.9865 0.941176 0.95679
## Neg Pred Value 1.0000 0.999698 0.99246
## Precision 0.9865 0.941176 0.95679
## Recall 1.0000 0.969697 0.86592
## F1 0.9932 0.955224 0.90909
## Prevalence 0.1306 0.009865 0.05351
## Detection Rate 0.1306 0.009567 0.04634
## Detection Prevalence 0.1324 0.010164 0.04843
## Balanced Accuracy 0.9990 0.984547 0.93186
## Class: Stroma Class: Treg Class: Tumor
## Sensitivity 1.00000 0.96667 0.9907
## Specificity 1.00000 0.99877 0.9930
## Pos Pred Value 1.00000 0.95604 0.9913
## Neg Pred Value 1.00000 0.99908 0.9924
## Precision 1.00000 0.95604 0.9913
## Recall 1.00000 0.96667 0.9907
## F1 1.00000 0.96133 0.9910
## Prevalence 0.03288 0.02691 0.4481
## Detection Rate 0.03288 0.02601 0.4439
## Detection Prevalence 0.03288 0.02720 0.4478
## Balanced Accuracy 1.00000 0.98272 0.9918Predict cell labels of unlabelled cells. Cells for which the highest class probability is below 40% are labeled as “unknown”.
cur_mat <- t(assay(unlab_spe, "exprs")[rowData(unlab_spe)$use_channel,])
cell_class <- as.character(predict.train(rffit,
newdata = cur_mat,
type = "raw"))
names(cell_class) <- rownames(cur_mat)
cell_prob <- predict.train(rffit,
newdata = cur_mat,
type = "prob")
cell_class[rowMax(as.matrix(cell_prob)) < 0.4] <- "undefined"
cell_labels <- spe$cell_labels
cell_labels[colnames(unlab_spe)] <- cell_class
spe$celltype <- cell_labels3.8.3 Cluster and cell type visualization
Visualize cell type and cluster labels on the UMAP embedding for qualitatively assessing cell phenotyping.
p1 <- dittoDimPlot(spe, var = "celltype",
reduction.use = "UMAP_mnnCorrected", size = 0.2,
do.label = TRUE) +
ggtitle("Cell types on UMAP, integrated cells")
p2 <- dittoDimPlot(spe, var = "nn_clusters",
reduction.use = "UMAP_mnnCorrected", size = 0.2,
do.label = TRUE) +
ggtitle("Clusters on UMAP, integrated cells")
p1 + p2
Visualize mean marker expression per cell type and per cluster as heatmaps.
library(scuttle)
library(viridis)## Loading required package: viridisLitecelltype_mean <- aggregateAcrossCells(as(spe, "SingleCellExperiment"),
ids = spe$celltype,
statistics = "mean",
use.assay.type = "exprs",
subset_row = rowData(spe)$use_channel)
dittoHeatmap(celltype_mean,
assay = "exprs", cluster_cols = TRUE,
scale = "none",
heatmap.colors = viridis(100),
annot.by = c("celltype","ncells"))
cluster_mean <- aggregateAcrossCells(as(spe, "SingleCellExperiment"),
ids = spe$nn_clusters,
statistics = "mean",
use.assay.type = "exprs",
subset_row = rowData(spe)$use_channel)
dittoHeatmap(cluster_mean,
assay = "exprs", cluster_cols = TRUE,
scale = "none",
heatmap.colors = viridis(100),
annot.by = c("nn_clusters","ncells"))
3.9 Spatial analysis
Perform spatial community analysis. This method groups cells solely based on their location in the tissue by using a previously constructed spatial cell graph. We perform community detection separately for tumor and non-tumor cells. An igraph object is constructed from the colPair slot storing the spatial cell graph constructed by steinbock, and the Louvain algorithm64 is used to detect spatial communities.
library(igraph)
# Spatial community detection - tumor
tumor_spe <- spe[,spe$celltype == "Tumor"]
gr <- graph_from_data_frame(as.data.frame(colPair(tumor_spe, "neighborhood")),
directed = FALSE,
vertices = data.frame(index = seq_len(ncol(tumor_spe))))
cl_comm <- cluster_louvain(gr)
comm_tumor <- paste0("Tumor_", membership(cl_comm))
comm_tumor[membership(cl_comm) %in% which(sizes(cl_comm) < 10)] <- NA
names(comm_tumor) <- colnames(tumor_spe)
# Spatial community detection - non-tumor
stroma_spe <- spe[,spe$celltype != "Tumor"]
gr <- graph_from_data_frame(as.data.frame(colPair(stroma_spe, "neighborhood")),
directed = FALSE,
vertices = data.frame(index = seq_len(ncol(stroma_spe))))
cl_comm <- cluster_louvain(gr)
comm_stroma <- paste0("Stroma_", membership(cl_comm))
comm_stroma[membership(cl_comm) %in% which(sizes(cl_comm) < 10)] <- NA
names(comm_stroma) <- colnames(stroma_spe)
comm <- c(comm_tumor, comm_stroma)
spe$spatial_community <- comm[colnames(spe)]
plotSpatial(spe[,spe$celltype == "Tumor"],
node_color_by = "spatial_community",
img_id = "sample_id",
node_size_fix = 0.5) +
theme(legend.position = "none") +
scale_color_manual(values = rev(colors()))
The step took 0.07 minutes.
Perform cellular neighborhood (CN) analysis. CNs are tissue regions with characteristic cell type composition and represent sites of unique local biological processes and interactions. We first construct a spatial 20-nearest neighbor graph based on the cells’ centroids, to include a larger number of neighbors compared to the spatial cell graph constructed earlier using steinbock. The aggregateNeighbors function of the imcRtools package then computes for each cell the cell type fraction among its 20-nearest neighbors. Cells are subsequently clustered using k-means clustering to detect CNs. Finally, these can be spatially visualized and the cell type fraction per CN can be computed.
library(pheatmap)
spe <- buildSpatialGraph(spe, img_id = "sample_id", type = "knn", k = 20)
spe <- aggregateNeighbors(spe, colPairName = "knn_interaction_graph",
aggregate_by = "metadata", count_by = "celltype")
set.seed(220705)
cn_1 <- kmeans(spe$aggregatedNeighbors, centers = 6)
spe$cn_celltypes <- as.factor(cn_1$cluster)
plotSpatial(spe,
node_color_by = "cn_celltypes",
img_id = "sample_id",
node_size_fix = 0.5) +
scale_color_brewer(palette = "Set3")
mat <- colData(spe) %>% as_tibble() %>%
group_by(cn_celltypes, celltype) %>%
summarize(count = n()) %>%
mutate(freq = count / sum(count)) %>%
pivot_wider(id_cols = cn_celltypes, names_from = celltype,
values_from = freq, values_fill = 0) %>%
ungroup() %>%
select(-cn_celltypes)## `summarise()` has grouped output by 'cn_celltypes'. You can override using the
## `.groups` argument.pheatmap(mat,
color = colorRampPalette(c("dark blue", "white", "dark red"))(100),
scale = "column")
The step took 0.18 minutes.
Perform spatial context (SC) analysis. SCs build up on the concept of CNs and are regions in which the local biological processes of CNs interact, and where specialized biological events take place. We construct a second k-nearest neighbor graph with larger k (k=40) to include cells across a length scale on which biological signals could be exchanged. For each cell, the aggregateNeighbours function computes the fraction of CNs among its 40 nearest neighbors. The detectSpatialContext function sorts the CN fractions from high to low, and the SC of each cell is assigned as the minimal combination of CNs that additively surpass a user-defined threshold (here threshold=0.9). After filtering the detected SCs, we can spatially visualize them and represent SC interactions as a graph.
spe <- buildSpatialGraph(spe, img_id = "sample_id",
type = "knn",
name = "knn_spatialcontext_graph",
k = 40)
spe <- aggregateNeighbors(spe,
colPairName = "knn_spatialcontext_graph",
aggregate_by = "metadata",
count_by = "cn_celltypes",
name = "aggregatedNeighborhood")
spe <- detectSpatialContext(spe,
entry = "aggregatedNeighborhood",
threshold = 0.90,
name = "spatial_context")
spe <- filterSpatialContext(spe,
entry = "spatial_context",
group_by = "patient_id",
group_threshold = 3,
cells_threshold = 100,
name = "spatial_context_filtered")
plotSpatial(spe,
node_color_by = "spatial_context_filtered",
img_id = "sample_id",
node_size_fix = 0.5,
colPairName = "knn_spatialcontext_graph")